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Abstract 

The paper begins with a short overview of the recent work done in the field of discontinuous reinforced 
composites, focusing on the different parameters which influence the material behavior of discontinuous 
reinforced composites, as well as the various analysis approaches undertaken. Based on this overview 
it became evident, that in order to investigate the enumerated effects in an efficient and comprehensive 
manner, an alternative approach to the computationally intensive finite-element based micrc^mechamcs 
approach is required. Therefore an investigation is conducted to demonstrate the utility of utilizing the 
generalized method of cells (GMC). a semi-anahtical micromechanics-based approach, to simulate the 
elastic and elastoplastic material behavior of aligned short fiber composites. The results are compa 
with 1) simulations using other micromechanical based mean field models and finite element (FE) unit 
cell modeLs found in the literature given elastic material behavior, as well as 2) finite element unit cell a 
a new semi-anahtical elastoplastic shear lag model in the inelastic range. GMC is shown to definitely 
have a window of applicability when simulating discontinuously reinforced composite material behav ior. 


1 Introduction 

Metal matrix composites (MMCs) have received considerable attention over the past 30 1 years due to 
their attractive specific strength, stiffness, fatigue and thermal properties[Clyne and V ithers (1993)]. [kelly 
and Z we ben (Eds) (2000)], [Sinclair and Giegson (1997)], [Lloyd (1994)]. MMC materials can be gener- 
icallv classified by their type of reinforcement, which consists of either continuous fiber (identified herein 
as CFMMC) or discontinuous particulate (termed herein as DRX) reinforcements. While both classesof 
materials share common features at a fundamental level (e.g„ matrix/reinforcement compatibility , rnterfa- 
cial properties sensitivity, architectural interaction (reinforcement size, shape and packing) and property 
opthnization. at a practical level dramatic differences can and do exist. For instance, CFMMC typically 
have exceptionally high strength and stiffness capabilities that are inherently strongly anisotropic; whereas, 
DRX material properties are ty pically more isotropic in nature with lower specific strength, y et have excel- 
lent stiffness capability . Clearly, failure to appreciate and consider explicitly the appropriate reinforcement 
morphology can result in inaccurate performance assessments, as elaborated upon m the next section. e 
combination of good transverse properties, relatively low cost, high workability and reasonable increases m 
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performance over monolithic metallic alloys have made discontinuous reinforced materials the most commer- 
cially attractive system to date for man) stiffness driven applications [Chile and Withers (1993)1. (Kelly and 
Zvveben (Eds) (2000)], [Miracle and Maruyama (2000)]. 

Although manufacturing issues and cost are probably the greatest impediments to the wide application of 
these advanced materials, it is also dear that accurate design and life prediction tools are greatly needed to 
facilitate the implementation of these developing materials. Although closure has not been reached regarding 
the best models for use in design and life prediction, it has become evident that if a model is ever to serve a 
purpose beyond that of basic research, it must fulfill several priman requirements. These indude a significant 
evel of accuracy on both the macro and micro scales, computational efficiency, and compatibility with the 
finite element method Fulfillment of these requirements allows a model to serve the composite de\eloper 
by enabling quick and easy variation of composite parameters for material development and optimization 
purposes. Likewise, those w ho attempt to utilize composites in structural design are well served if the model 
is compatible with the finite element method. The generalized method of cells (GMC). originally developed 
b ; Aboudi (1991. 1995), is a good choice for implementation into modeling tools for advanced composites, 
requirem ents described above. Particularly, whenever microfield as well as macrofield quantities 
w ill be needed - as is the case when attempting to account for inelastic material behavior in MMCs. 

GMC is a full) analytical micromechanics model for multi-phased materials with arbitrary periodic 
microstructures. It provides pseudo closed-form multiaxial constitutive equations for such materials, and 
allows straightforward implementation of physically-based visooplastic deformation models, as well as arbi- 
trary failure and damage models for each phase. Further, recent independent advances have simplified the 
implementation of GMC as an elemental constituent material model in finite element analysis [Arnold et al. 
(1999)], and significant!) increased the model's computational efficiency [Pindera and Bednarcv k (1999)]. 
CMC has been implemented in the NASA Glenn Research Center comprehensive micromechanics analvsis 
code, MAC/GMC [Arnold et al, (1999)]. The code has man) features that render it useful for design, de- 
formation modeling, and life prediction for a wide range of materials. These features include the ability to 
simulate general thermomechanical loading on composites whose geometries are represented bv a library of 
continuous and discontinuous repeating unit cells, a library of nonisothermal elastic/ viscoplastic constitutive 
models, fatigue damage anal) sis, yield surface anal) sis, laminate analysis, and interface modeling. 

In the past, the various capabilities described above have been exercised extensively, however, the pre- 
dominant material system examined has been that of the CFMMC [Aboudi(1996)], [Arnold, et. al. (1996)1. 
[Goldberg and Arnold (2000)]. [lyeretal. (2000)], [Bednarcv k and Arnold (2000a)], [Lissenden et al. (2000)] 
[Bednarqyk and Arnold (2000b)]. Consequently, the priman objective of this current stud) Is to assess the 
applicability of GMC in predicting the elastic and inelastic (plastic) behavior of particulate reinforced com- 
posites (DRXs). The literature indicates that DRX behavior is significantly influenced bv various charac- 
teristic parameters like fiber arrangements (staggered, non^taggered, cluster, random, packing type), fiber 
volume fraction, fiber aspect ratios, fiber-matrix stiffness ratios and reinforcement type (cylinder, particle, 
etc. ). A brief review- of this is provided in the background section. Even though GMC is full) capa- 
ble of modeling arbitrary multiphased composite architectures, to facilitate comparison w ith previous work- 
reported m the literature and limit the scope of this investigation, only one type of DRX is considered, 
namely, the simulation of aligned short fiber composites with non-staggered cylindrical reinforcements and a 
hexagonal packing arrangement. The need for such an assessment stems from the fact that although GMC 
consists of analytical equations that are able to deal with both the homogenization and localization process 
for composite materials in a very numerically efficient way, the method’s linear displacement field assumption 
in combination with the imposition of displacement and traction boundary conditions in an average sense, 
gives rise to an inherent lack of normal and shear component coupling. This lack of coupling is such that 
application of a global normal (or shear) stress field will produce only a local normal (or shear) stress field. 
Consequently, for cases (e.g.. DRX materials) in which the load transfer mechanism between phases (e.g.. 
matrix and fibers) is shear dominated, the lack of normal-shear coupling may become quite problematic 
In particular, in the case of an elastoplastic matrix regime placed between two short fibers, GMC’s lack of 
shear coupling requires both the fiber and matrix to cam the same normal stress component. This tends 
to produce an overly compliant (soft) composite response. The results presented within demonstrate this 
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fact and more importantly characterize the window of applicability for the GMC method within the context 
of discontinuous!)' reinforced metallic composites. The introduction of an "artificial-’ elastoplastic matrix 
material between the particulate reinforcements. Ls proposed in an attempt to extend the applicability of 
GMC. Although, providing more reasonable macro behavior, such a strategy still results in incorrect mi- 
cromechanical stress fields within the unit cells, thereby limiting the usefulness of the current GMC model 
when calculating failure or fatigue analy ses. 

With section 2, the paper begins by reviewing recent work done in the field of short fiber composi es so 
as to reveal the kev parameters influencing DRX material behavior. Section 3 then articulates the various 
analytical and numerical approaches employ ed and compared in this paper. Finally , section 4 presents a 
comparison of GMC results with other analytical and numerical results within the elastic and elastoplastic 

range, respectively. 


2 Background: Factors Influencing DRX Materials 

The early work on the principles of strengthening in composites reinforced with discontinuous fibers, 
whiskers or particulates carried out in the 1950's and 1960 s Ls characterized by the use of relatively simple 
models that do not explicitly take into account the actual inclusion distribution through interaction between 
adjacent inclusious. Many modeLs have been employ ed to study the effect of inclusion shape on the elastic 
moduli, as well as subsequent inelastic response, and also the mechanism of stress transfer from the matrix 
to a finite-length fiber. For instance. Eshelby’s analytical .solution to the problem of an arbitrarily shaped 
ellipsoidal inclusion embedded in an homogeneous material deformed by uniform tractions or displacements 
at infinity has formed the basis for calculating the effective response of macroscopically homogeneous two- 
phase composites using a number of approaches [Eshelby (1957)]. These approaches include the self-cousistent 
scheme which neglects the inclusion/ matrix interaction in calculating stress fields in the inclusion phase [Hill 
(1965)], and the Mori-Tanaka method (MTM) which takes this interaction into account in an approximate 
fashion [Mori and Tanaka (1973)]. The problem of an array of ellipsoidal inclusious with different aspect 
ratios embedded in an elastoplastic matrix can be treated using the Mori-Tanaka approach, as was done 
bv [Brown and Clarke (1975)] investigating the effect of inclusion shape on work hardening of metal matrix 
composites. However, this method ty pically underestimates y ielding and subsequent hardening effects due to 
the use of mean stress and strain fields in the matrix phase, and treats all inclusion distribution on the same 
footing so long as macroscopic homogeneity is preserved. Along similar lines, the so-called shear-lag analysis 
has been employed to study the effectiveness of short fibers as reinforcement using the strength-of-matenals 
approach to analy ze stress" fields around and within a finite-length fiber embedded in a surrounding matrix 
[Dow (1963)]. While this type of analy sis helps to identify shearing of the matrix as the primary mechanism 
of force transfer from the matrix to the fiber, and thus the critical fiber length over which the axial stress is 
introduced into the fiber from both ends, it Ls based on a very simplified one-dimensional analy sis of stress 
fields which neglects the influence of morphology of surrounding fibers, among other things. 

Despite the relatively long history of modeling the response of discontinuous fiber composites, only 
recently have sy stematic investigations of the effects of reinforcement shape and arrangement been initiated 
for these tv pes of composites. Inclusion of the third dimension in the analy sis of the elastic and inelastic 
response of discontinuous fiber composites increases the number of variables several fold relative to the 
analy sis of continuous fiber composites. Thus, in addition to the arrangement and shape of the reinforcement 
in the plane transverse to the loading direction, variables associated with the planes parallel to the loading 
direction were included. These variables include the fiber aspect ratio, fiber spacing in the vertical and 
horizontal directions, those associated with the effect of fiber clustering, and the extent of overlap between 
adjacent columns of fibers. These added complexities typically require numerical solution procedures, such as 
the finite-element analysis: particularly if complicated reinforcement shapes and arrangements are mvolved. 

Finite-element investigations carried out in recent years have focused on separating the effects of inclusion 
shape from the effects of inclusion packing array geometry on the overall elastic, elastoplastic and creep 
behavior. Due to many different tv pes of reinforcement employed in discontinuously-reinforced composites, 
including whiskers, particulates, platelets as well as finite-length fibers, extensive analyses of inclusion shape 
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have been conducted for this class of composites. Examples of inclusion shapes investigated include spherical 
particles, circular cylinders with different aspect ratios, truncated circular (octagonal) cylinders, double-cone 
particles, and ellipsoidal or disk-like particles (platelets). Similarly, diflerent packing arrays studied include 
hexagonal arrangements of inclusions in the plane of loading distributed in either regular (aligned) arrays or 
off-set (staggered) arrays in the planes parallel and perpendicular to the loading direction so as to model and 
assess the influence of distribution randomness that is typical in DRX composites. Inclusions such as disk- 
like particles arranged in packet morphologies have also been considered, as have clustered arrangements. 
These investigations have been driven in large part by the wide range of inclusion morphologies that result 

from current material processing techniques, as well as by the wide range of shapes available for the inclusion 
phase. 

In order to reduce the complexity of a three-dimensional finite-element formulation in analyzing the 
response of discontinuous fiber, w hisker or particulate composites, various idealizations of the unit cell have 
been employed by a number of researchers (e.g. clusters of whiskers or cylindrical particles in the plane 
parallel to the applied load have been modeled using a plane strain idealization of rectangular platelets). 
This effectively reduces the problem to that of a continuously reinforced composite subjected to loading 
in the plane perpendicular to the long fiber direction. The results of such analysis for rectangular cross- 
sections with different aspect ratios in the plane of loading .should be easily deduced from the early analysis 
of continuously reinforced composites with rectangularly shaped fibers [Ashton et al.(1969) ]. In the case of 
a hexagonal array of inclusions in the plane transverse to the loading direction, the problem is often reduced 
to an axisymmetric (i.e., two-dimensional) problem by approximating the unit cell using a circular cy linder 
with different types of lateral boundary conditions to simulate the interaction with adjacent fibers (i.e., unit 
oells). Usmg this model, different fiber arrangements in the plane parallel to the applied load have been 
investigated, with varying amounts of overlap between vertical columns of adjacent fibers (e.g. [Christman 
et al. (1989)], [Tvergaard (1990)], [Dragone and Nix (I960)], [Povirk et al. (1990)], [Yang et al. (1991)] 
[Bao et al.(1991)], [Siegmundet al. (1992)]). In contrast to continuously reinforced composites, the inclusion 
arrangement in discontinuous^ reinforced composites is inherently three-dimensional, so that the use of 
plane models requires considerable caution. The results of fully three-dimensional finite-element analyses 
of such periodic array s are also available (e.g. [Levy and Papazian (1990)], [Weissenbek and Rammeretorfer 
(1993)], [Abel et al.(1993)], [Weissenbek et al. (1993)]). They are, however rather restrictive in terms of 
fiber arrangements and shapes, and tend to entail considerable computational costs. Both three-dimensional 
and axisymmetric unit cell models have been used successfully for studying the nonlinear thermomechanical 
behavior of aligned short reinforced composites. 

The literature survey presented indicates significant influence of both inclusion shape and distribution on 
the elastic, elastoplastic and creep response of discontinuous fiber composites. The extent of this influence for 
a given fiber architecture depends on the inclusion content, inclusion/matrix material property mismatch, and 
the direction of applied load with respect to the internal micro-structure. The influence of fiber architecture 
on the response of DRX oomposites is significant and is based on a substantial number of investigations, 
each of which was necessarily limited in scope due to the prevalent use of the finite-element approach in 
modeling the response of the investigated architectures. Thus while a considerable body of knowledge 
has been generated that sheds fight on the effectiveness of different fiber architectures in strengthening 
discontinuous fiber composites, considerably more sy stematic research Is required to develop design guidelines 
for optimization of material performance through fiber architecture 

Hovvever. it becomes dear, based on the above discussion, that in order to investigate the enumerated 
effects in an effident and comprehensive manner, induding parametric studies involving fiber content and 
material property mismatch variations, an alternative approach to the computationally intensive finite- 
element based micromechanics approach Is required This is particularly true when inelastic (viscoplastic) 
theories, which typically require computationally intensive integration algorithms, are employed to model 
the response of metallic matrices. Therefore, the objective of this paper is to apply the computationally 
effident GMC method so as to determine its suitability for predicting the macroscopic response of aligned 
DRX materials. 
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3 Modeling Approaches 

All simulations conducted in this study consist of aligned short fiber composites, wherein the particulates 
(wishers) am idealized as non-staggerod. cylindrical reinforcements, with a hexagonal packing arrangement. 
S imilar ly other basic constraints imposed are as follows: 


The inclusion properties are assumed isotropic and linear elastic, while the matrix is taken to be 
isotropic elastoplastic with isotropic hardening. The material data for both the inclusion and matnx 
are treated as temperature-independent and isotropic. 


• The inclusions are axis: mine trie, identical in shape and size and can be characterized by an aspect 
ratio, « r . that relate the length of an inclusion to its diameter. 

• Investigated model parameters are the fiber (inclusion) volume fraction (V f = 0.2. 0.5), the fiber-matrix 
stiffness ratio (s r = E F /E M = 3. 30) and the fiber aspect ratio (a r = l/cl = 5. 15. -5). 

• The inclusions and matrix are well bonded at their interface and remain that way during deformation. 
Thus, we do not consider interfacial slip, fiber-matrix debonding or matnx nncro^rackmg. 

• No time dependence (relaxation or creep) effects are included in the present analysis. 


Note, obtaining experimental reference data for unidirectional short fiber composites appears to be problem- 
atic. as it has not proved possible to produce physical samples with perfectly aligned fibers. Consequent 1 } , m 
this paper we avoid this additional complexity of orientation effect by using threeehmensronal finite element 
analysis of aligned short fiber oomposites as our reference standard, rather than experimental results. 


3.1 Material Data 

The material properties used for all calculations are given in Table 1. These properties are in the range 
of typical fiber-reinforoed engineering thermo-plastics (PMC’s - Material 1) and metal matnx composes 
(MMC’s - Material 2), respectively. The following input data are needed for both components^ oungs 
modulus E and Poisson s ratio v. The yield stress a Y and the hardening modulus E H are needed for the 

matrix material only. 


Table 1: Prototypical material properties used for all calculations. 



Material 1 (sr =30) 

Material 2 (s^ =3) 

Property 

Fiber 

Matrix 

Fiber 

Matnx 

E (GPa) 

300. 

10. 

300. 

100. 

V 

0.17 

0.33 

0.17 

0.33 

cry (MPa) 

- 

20. 

- 

200. 

E h (GPa) 

- 

0.1 


1.0 


3*2 Analytical Models 

In the past, numerous analytical models have been examined in an attempt to determine the ‘best ’ 
(see e.g., [Tucker and Liang (1999)], [Aboudi (1991)]) available analytical model for predretmg the elastic 
stiflnefs behavior of aligned short fiber reinforced composites relative to reference FEA-mode 
tations. Analytical models considered, were modeLs such as: the Mon-Tanaka method (MTM)[Mon and 
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Figure 1: Sketch of inclusions-matrix shapes for the applied aligned short fiber composite models, (a) Mori- 
Tanaka model, (b) Hashin-Shtrikman estimates of Ponte Castaneda and Willis, (c) GMC-model (a r = 5.0) 


Tanaka (1973)], self-consistent modeLs. bounding models (e.g., variable constraint model (VCM) and Hashin- 
Shtrikman bounds and estimates), Halpin-Tsai equations and shear lag modeLs (see Appendix B for details 
of the elasto plastic extension used in this study). For example, the recent work of [Tucker and Liang ( 1999)] 
recommended that the Mori-Tanaka model was the “best” choice for estimating the overall elastic stiffness 
behavior of aligned short fiber composites. 

In this stud\- we will extend the scope of comparison to include the semi-analrtical generalized method of 

^rT < ^ I rl GIVIC) a,ld !° mpare GMC " ith the WTM, VCM, FEM results and Hashin-Shtrikman estimates 
(HSE) of [Ponte Castaneda and Willis (1995)]. Benveniste (1987) has provided a particularly simple and 
clear explanation of the Mori-Tanaka approach. Analytical bounds for the stiffness behavior can be obtained 
from the variable constraint model (VCM). published by [Pedereen and Wither (1992)]. Freni a practical 
point of view it is of interest that the lower (VCM) bound corresponds to the standard Mori-Tanaka results, 
while the upper (VCM) bound can also be obtained from the Mori-Tanaka approach after a so called “color 
inversion” (i.e.. exchanging the roles of inclusion and matrix). These bounds obtained from the Mori-Tanaka 
4130 correspond to the Willis bounds in the case of aligned inclusions, compare for example, [Weng 

It is important to remember that the basic assumption of the Mori-Tanaka method is that the aligned 
ellipsoidal inclusions are surrounded by coaxial aligned ellipsoidal matrix domains in which both ellipsoids 
have the same aspect ratio (Fig. la) and that each inclusion (fiber) within a concentrated composite sees the 
average strain of the matrix. Consequent!} , the MTM provides only macro (average) fields and properties 
and therefore lacks the ability to address adequately problems dominated bv microfield quantities, such as 
damage and inelastic behavior anal} sis. Alternatively, GMC allows both homogenization and localization of 
the associated stress and strain fields (see Appendix A), thereby providing not only macro stiffness properties 
but also local eigenstrain field quantities. Note, that in the present investigation the GMC unit cell differs 
*4™ 88 Afferent fiber-unit cell aspect ratios (Figs.2 and 3) can be realized since the whole inclusion 
(fiber) is surrounded by a matrix material layer of constant thickness (side-to-side and end-to-end distances 
are equal) and labeled, e. subsequently. Because of this varying aspect ratio, HSE for the prediction of the 
effective moduli of the elastic inhomogeneous material given in [Ponte Castaneda and Willis (1995)] are also 
introduced (Fig. 1 b) to overcome the drawback of the original Mori-Tanaka method so that different fiber-unit 
cell aspect ratios can be realized Note, these Hashin-Shtrikman estimates correspond identically with the 
Mon-Tanaka method for the case of aligned coaxial inclusions with equal aspect ratios of the inclusion and 
the surrounding matrix ellipsoids (see [Hu and Weng (2000)]) 


J €3 ‘ imateS differ t>om the Hashin-Shtrikman bounds, as [Ponte Castaneda and Willis (1995)1 

[ lllis (19|7)] removed the restriction of statistically isotropic overall material symmetry (used bv both HS and Willis) 

^r rat T * j* Sp , at,al distribu ‘ ion of ‘"elusions from the inclusion shape. Also in the HSE. [Ponte' Castaneda and Willis 

17k S 3 SUl ® reterence raaterial th * effective material like in self consistent models) instead of the two 

used by [Willis (i9f/)j. 
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3.3 Generalized Method of Cells 

Although, GMC Is capable of analyzing any multiphased composite material with an arbitral? internal 

microstructure and reinforcement shape (see Appendix A), the rw m^te 

in Fig 2. Due to the periodic boundary conditions inherent to the GMC-model the whole R\ E must be 
modeled, instead of only a quarter of the unit cell (see the dashed rectangle in Ft* 2) as typically done 
in the case of FEA representations where symmetric boundan conditions are applied. Apart from this 



Figure 2: Aligned short fiber GMC unit cell model, comprised of 338 subcells. 

fact the unit cell dimensions were chosen similar to the dimensions of the FE-model (see Fig. 3). Th «J fib * r 
was ' approximated with more than one subcell (sufficient for the longitudinal stiffness behavtor) m order to 
reduce the errer in the transverse response (i.e., the difference between £22 and £ 33 ) and better approximate 
the targeted transvereeh isotropic response behavior resulting from a hexagonal packing arrangement. The 
cross sectional area of the approximated fiber Ls taken equal to the cross sectional area of the actual circular 
fiber. These conditions, together with the known fiber diameter, fiber volume fraction, and the distance, e 
(see Fig. 3) suffice to determine the dimensions of the unit cell for each fiber 

GMC analyses were realized with the computer code MAC/GMC [Arnold et al.(1999)] veoped a e 
NASA Glenn Research Center. 

3.4 Finite Element Models 

The finite element unit cell analyzed consisted of an identical periodic, three dimensional non-^taggered 
array of fibens (similar to that idealized with GMC) with the corresponding dimensions ffiown in Fi«- 3- The 
unit' cell dimension were chosen such that b would be proportional to ft (i.e., b = Jo) through the constant 
3. In the case of hexagonal packing, 3 = V5, and the distance between neighboring fiberends ^-to W 
distance, equal to c - l f in Fig. 3) was set to the closest distance between two fibers (side-to-side distance, 

equal to e in Fig. 3). 

These conditions, together with the fiber diameter and volume fraction, suffice to detennine the dimension 
a.b and c for the repeating volume element (RVE) to be analyzed. Note that a new unit cell and its 
corresponding 3D mesh are generated for each fiber aspect ratio (l/d) and fiber volume fraction {v f ). 
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GMC unit cell 


FE unit cell 




Figure 3: Example representative volume element and finite element unit cell for a h = 5.0 and Vp = 0.2 

Stiffnesses within the elastic and plastic range of these unit cells were calculated using [ABAOUS (199S)1. 
ully mtegrated. first order. soUd elements 2 were used, as this element type does not sufler from volumetric 
locking, which is critical for accurate elastoplastic calculations. The mesh discretization containing 13.300 
elements employed is shown in Fig. 3. This mesh density is approximately three times finer than that used 
; ^ e f Sen i bek fvanunerstorfer (1993)] and was shown during the course of this work (by convergencv 
studtes) to be sufficient!} accurate for both elastic and inelastic analysis. Modeling of only one eighth of the 
w ole R\ E (depicted on the left) was achieved due to the use of appropriate symmetric boundary conditions. 

or example, when considering axial or transverse loading, symmetry requires all faces of the unit cell to 
remain plane. Consequently, to determine E u and u\ 2 the normal displacements of the front (xj = 0). left 
(x 2 - a, -) and bottom (x 3 = 0) faces of the unit cell were fixed Then a master node was defined (the 
intersection point of the back, right and top face) and linear constraint equations applied on the back, right 
and top face such that the normal displacements of all nodes (except the master node) on these faces were 
set equal to the normal displacements of the master node. The tangential displacements on all faces were left 
unstrained The average (marco) stress was then computed from the reaction force of the master node 
vided by the cross-sectional area normal to the reaction force of the investigated unit cell. Similarly the 
average (macro) strain w as computed by dividing the displacement of the master node bv the corresponding 
length of the unit cell. For convenience, a one Newton load was applied to the unit cell 

Analogous conditions were used to load the unit cell in the Erection to determine E 2 , and 
The fifth and final independent material parameter for transversal isotropic material, the longitudinal shear 
moduhis C n2 . could also be determined using the unit cell shown in Fig. 3, assuming the appropriate periodic 

boundary conditions were applied to each face. This more complicated analy sis was not undertaken at this 
time. 


4 Results 

4.1 Elastic Material Behavior 

AU of the analy tical micromechanics theories used in this study predict full transversely isotropic prop- 
erties. Transverse isotropy about the xj ax is implies that the stillness behavior is the .same for any loading 
2 In [ABAQUS (199S)] this is the element known as C3D8. 
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direction in the 2-3 plane; consequently, this not only requires that E 22 = E zz . but also that 


£*23 — 


2(1 + "23 ) 


Therefore, onlv fire independent stiffness constants (e.g. , E n . E-a, v 12 , "23 and G 12 ) need be determined or 
the case of transverse isotropy. Periodic unit cells with hexagonal packing should be transversely tropic 
as there are six planes of symmetry. However, since periodic boundary conditions were not applied to e 
FE unit cell, we were unable to obtain Gy 2 from the FE investigations. Consequently, we will limit our 
comparison for the varions approaches to the prediction of Eu, E 22 . "12 and y 2Z only. 

Results of the normalized axial modulus (xi - direction, .scaled by the matrix modulus) obtained from e 
VCM-model (MTM-model), Hashin-Shtrikman estimates, GMC-. and FE-model are shown in ig. ( 
like materials with a lower stiffness ratio, e.g.. s r = 3) and Fig. 5 (PMC like material, with a higher strips 
ratio eg s r = 30) for numerous fiber aspect ratios and two fiber volume fractions, a low (v/ - 0~) and 
a hi A \v f = 0.5). The Hashin-Shtrikman-, CMC- and FE-modeLs wore evaluated only at discrete aspect 
ratios; whereas, the VCM-result were computed for many different aspect ratios and thus displayed wi a 
solid line. Note, for high stiffness ratios only the lower VCM-bound (which corresponds to the hITM) was 
plotted. From Figs. 4 and 5 it is apparent that the influence of fiber aspect ratio on Eu is sigmfican ( ax. 
error of 35%) for a composite with a high stiffness ratio whereas relatively insignificant (Max. error of 0 / 0 ) 
for materials with a low stiffness ratio. In both cases the Hashin-Shtrikman estimate lay asexpected vvithm 
the upper and lower VCM bounds. The GMC- and FE-results are partially out of the Mon-Tanka based 
VCM-bounds due to the different inclusion topology (ellipsoidal versus cylindrical). Generally spiking, 
the GMC-model underestimates the longitudinal stiffness, particularly for high stiffoess ratios, although the 
trend is captured quite nicely. Alternatively , for materials with low stiffoess ratios, good agreement between 
GMC and Hashin-Shtrikman estimates as compared with. FE-results are observed given a \\i range o 
aspect ratios and fiber \olunie fractions. 


— VCM L.B. 

— VCMU.B. 

* HS-Estimates 
O FE-Results 
A GMC-Results 


V =0.5 

F 


v =0.2 

F 


Aspect Ratio, 


Figure 4: Comparison of axial modulus ratios for s r — 3 

The physical reason for the difference between the FEA-results and those of GMC can be inferred from 
Fig. 6, where the longitudinal stress distributions for a load in the x, -direction are shown. In both FEA 
analyses shown in Fig. 6, the geometry, boundary conditions, and applied forces are the same. Only e 
stiffness ratio between matrix and fiber constituent material is different. Evidently , these stress contours 
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Figure 5: Comparison of axial modulus ratios for s r = = 30 

indicate that an increase in the stiffness ratio leads to a decrease in the matrix loading and conversely to 
an increase in the fiber loading (see maximum and minimum stress values). Similarly , the difference of the 
"mean matrix stress in the volume of material between the two fiber ends as compared with the ‘•mean” fiber 
stress becomes greater as the stiffness ratio Ls increased. Considering, the GMC representation, the stresses 
are constant and equal along a longitudinal fiber-matrix row. Therefore the error for lower stiffness ratios 
is decreased, yet when the stiffness ratio increases so too does the error. Figure 6 also gives an explanation 
as too vvhy the CMC stiffness results are slightly under the FE-results in Fig. 4. as the GMC assumptions 
lead to an unloading of the “stiff” fiber and to an additional loading of the "soft” matrix. 

Figures 7 and 8 illustrate the normalized transverse modulus results obtained from the various analytical 
and numerical methods considered. Here as one would expect, the fiber aspect ratio has little, if any. 
effect on the transverse modulus Also, in general the GMC results compare favorably (maximum error of 
approximately 5%) with the FE-results throughout. For both stiffness ratios, the Hashin-Shtrikman estimates 
are equal to the lower VCM bound (which in turn are equal to the original MTM-results). Noteworthy is the 
fact that the GMC-results are much closer to the FE-results for the case of high stiffness ratios, low" aspect 
ratios, mid high fiber volume fractions than the other analy tical estimates. Similarly, as observed for the 
longitudinal stiffness the qualitative trend of the FE results is better captured by GMC. 

To help explain this observation, we examine the calculated transverse stress distribution obtained from 
FEA, as shown in Fig. 9. Note, only the case for the low’ stiffness ratio (sr = 3) Ls shown in Fig. 9 as 
the FE-results for the higher stiffness ratios appear similar. Relative to the FE-results for the longitudinal 
loading case, the difference in the "mean” matrix stress of the volume of material between the two fibers 
m the loading direction and that of the "mean” fiber stresses for the transverse loading case are small. 
Consequently , the inaccuracy caused by having similar stress states in a given row’ or column of subcells 
within GMC is small; thus explaining the lack of influence the stiffness ratio has on the GMC predictions 
relative to those determined using an FEA unit cell. 

Figures 10 and 1 1 show a comparison of the axial Poisson’s ratios based on the various evaluation methods. 
Apun the comparison between FEA and GMC results are quite favorable, irrespective of aspect ratio and 
stiffness ratio: the difference being approximately a maximum of 5.2%. Again, the GMC results follow the 
qualitative trend of the FE-results much better than either the Hash in Shtrikman estimates or MTM-results, 
especially for the case of high stiffness ratios and low’ fiber aspect ratios. 
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Figure (i: Example for the longitudinal stresses obtained from FE-calculation for s r - 3 and 
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Figure 9: Example of transverse stresses obtained from F E-calculation for s r = 3 (V F = 0.2, an - 5.0) 
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Figure 10: Comparison of axial Poisson ratios for s r = 3 
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Figure 11: Comparison of axial Poisson ratios for s r = 30 

Similarly, results for the in-plane Poisson’s ratio i/ 23 appear in Figs. 12 and 13. Once again, in general 
the GMC results are slightly above those produced with FEA. The Hashin-Shtrikman estimates being closer 
to the FE-results than the GMC results. The discrepancy between methods being more evident with higher 
fiber volume fraction than higher stiffness ratios. Even so, the maximum error between the GMC predictions 
and those of FEA is at most 4.8%. 

4.2 Elastoplastic Material Behavior 

In this section the elastoplastic stiffness behavior of aligned short- fiber composites will be analyzed using 
primarily the GMC and the finite element method. A comparison of the resulting predictions should confirm 
whether or not the GMC is capable of accurately simulating the elastoplastic stiffness behavior of DRX 
materials. Additionally, results obtained from a semi-analytical elastoplastic shear lag model (see Appendix 
B) are included as well to help illuminate any differences between the methods. The investigated model 
parameters are similar to those of the elastic investigation, with the exception that only two fiber aspect 
latios (ci r o and 25) are examined. Furthermore we restrict ourselves to the case of longitudinal (xj — 
direction) loading so as to allow the inclusion of shear lag results. 

Figure 14 shows the effective longitudinal stress-strain curve produced by the GMC (dotted and dashed 
line), FE unit cell (solid line), and shear lag (symbols) models for the special case of Vj = 0.2, a r = 5, ,s> = 3. 
Clearly, a large difference in the longitudinal stiffness response exists, depending upon the analysis method 
employed, with the GMC model predictions being significantly softer. The question before us is whether or 
not the lack of shear coupling inherent in the GMC method is the primary cause for the softer elastoplastic 
behavior. To address this question, we introduced a semi-analytical shear-lag-model (and shear-free model) as 
described in Appendix B. A comparison of the shear-lag (open circles) and shear- free (open square symbols) 
models with those of FE and GMC are also shown in Fig. 14. 

Evidently, the 8-subcell (1 fiber subcell surrounded by 7 matrix subcells) GMC results and the shear-free 
model results agree extremely well; whereas, those determined using the elastoplastic shear-lag model agree 
with those coming from the FE unit cell model. Consequently, it may be concluded that: 

1. The difference between the GMC- and FE-simulations is directly related to the lack of normal-shear 
coupling inherent in the GMC formulation. 
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Figure 12: Comparison of transverse Poisson ratios for s r = 3 



Aspect Ratio, a R 

Figure 13: Comparison of transverse Poisson ratios for s r — 30 


NASA/TM — 2001-211 165 


15 






Figure 14: Comparison of FE-. CMC-, Shear-lag-, Shear-free-model results (s r = 3, Vj- = 0.2, a r = 5) 


2. The standard (without any modifications) GMC-model is not well suited for predicting the elastoplastic 
response behavior of DRX materials, as the GMC-model predicts an overl.v .soft plastic response. 

3. The semi-analytica! elastoplastic shear-free and shear-lag models introduced in Appendix B appear 
folly capable of predicting both the GMC- and FEA-results. The primary disadvantage of this shear- 
lag model being its uniaxial nature, that Is, it can only provide the effective modulus (response) in 
a single direction. 

In addition to the pseudo-square 3 fiber (8-subcell) GMC representation (dotted line) in Fig. 14, the 
typiad (338 subcell) circular fiber unit cell (dashed line) Ls shown in Fig. 14. Clearly , a discretization effect 
is observed and is a result of the piecewise purification of the matrix subcells between the two fiber ends 
m the case of the finer GMC-model. This effect results in a softer response within the initial plastic range 
(i.e., from 0.2 - 0.5% strain range). After purification of all matrix subcells the effective hardening slope of 
both GMC unit cells are essential!} the same. 

To verily the above conclusions, simulations with other aspect ratios, volume fractions and stiffness ratios 
were conducted with the results being displayed in Figs. 15 - 18. The semi-analytical shear-free model is not 
shown m these figures as its results alway s agreed with the GMC (pseudo square) simulations. The fitting 
parameter v of the semi-analy tical shear-lag-model is found to be only a function of the fiber volume fraction 
i.e. v = w(Vf) and takes a value of approximately 2 for V) = 0.2 and approximately 2.5 for V f = 0.5. The 
dashed lines in all figures represent the effective stress-strain response of the standard circular GMC-model. 
All figures confirm the abo\e conclusions. 


Some characteristic features of the effective finite element stress-strain curves in Figs. 15 - 18 are evident- 
for example, increasing either the aspect ratio or volume fraction significantly hardens the composite response, 
particularly for high stiffness ratio materials. The elastoplastic stiffness response may be partitioned into 
three ranges. Figure 19 shows these ranges as a function of effective longitudinal strain. 

Range 1, delineates the domain where linear elastic behavior is manifested This Ls followed by range 2 
m which a large variation m the hardening slope takes place and range 3 in which approximately a constant 

*ln GMC the geometry of the subceU representing the fiber is square: however in actuality since all fields are taken at the 
centroid of the subcell only the ratio of subcell to total length is important, no comers are felt. Thus the square fiber is reallv 
a pseudo-square. 
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Figure 15: Comparison of FE-, GMC- and Shear-lag-model results for s r — 3 and 1 / — 0.2 



Figure 16: Comparison of FE-, GMC- and Shear-lag-model results s r = 3 and V/ = 0.5 
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Effective Longitudinal Strain, £ eff 


Figure 17: Comparison of FE-, GMC- and Shear-lag-model results s r =30 and Vf 



Figure IS: Comparison of FE-, GMC- and Shear-lag-model results s r = 30 and Vf 
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Figure 19: Semi-analytical shear-lag model, characteristic ranges and number of yielded truss elements 

(s r = 3, V f = 0.2, a r = 5) 

effective hardening slope Ls observed. More detailed investigations into the flow behavior show that at the 
beginning of range 2 the matrix material between the two fiber ends begins to yield first (see Fig. 19, point 
1). After which the whole matrix vields (point 2). During the longitudinal loading step the matrix from the 
side-to-side region tries to move to the end-to-end region. This material flow’ leads to a considerable shear 
leading of the side-to-side matrix interphase (between fiber and matrix). In the plastic range large shear 
deformations of the interphase start at the fiber ends (point 3), then grow’ over the fiber and prevent any 
additional loading of the fiber. In the case of the semi-analytical shear-lag-model this means that one bp 
one each non-linear truss element will start to yield All trusses are plastic at point 4 (see Fig. 19), which is 
the beginning of range 3. Note, the stair-stepped dashed line in Fig. 19 indicates the piecew ise yielding of. 
altogether, five truss elements. 

Based on these results, the question remaining is whether the GMC can be appropriate!} modified in 
order to predict accurately the elasto plastic behavior of DRX materials. Considering wh} the GMC fails, 
it becomes clear that if the matrix between the two fiber ends 3 ’ields, the fiber Ls no longer capable of 
earn ing anv additional load beyond the matrix since both the matrix and fiber have the same stress state. 
This continuity of specific stress components along a given row’ or column of subcells is a direct result of the 
inherent lack of normal and shear coupling with GMC. Given this fact, a simple modification to GMC would 
be the introduction of an artificial matrix material for all suboells within the region between the two fiber 
ends. This new’ artificial matrix material must be chosen so that the resulting effective stress-strain curve 
produced corresponds with the reference data (be it experimental or analytical). Of course this strategy 
the draw back that the obtained micro-mechanical GMC stress state is no longer based on anv 
physical meaning: i.e., fracture, damage or fatigue analyses could no longer be performed confidenth when 
using the modified GMC-results. as the microfield quantities would be even more questionable than before. 

As a first attempt, the artificial matrix material was assumed to behave as purely elastic, with the 
associated Young s modulus of the artificial matrix material being that of the matrix itself. Clearly, such an 
assumption results in the overestimation (at least well within the plastic range) of the elastoplastic response 
curve as illustrated in Fig. 20. This response prediction is non-conservative. In contrast, it Ls shown that the 


NASA7TM — 200 1-211165 


19 




Figure 20: Comparison of FE-, GMC-, and GMC-modified simulations (s r — 3. Vf — 0.2, a r = 5). 


standard GMC-model provides a lower bound (conservative response prediction) for the effective elasto plastic 
behavior, thus suggesting that an artificial bilinear elastoplastic material might provide a more accurate 
simulation. For such a material definition, four material parameters would be required as input; the} are 
1) Poisson s ratio, 2) Youngs modulus, 3) yield stress, and 4) the plastic hardening slope. With the yield 
stress and hardening slope being the best candidates for modification from that of the matrix material, since 
the elastic range is already accurately predicted, the parameters (E,/y) can be taken to be that of the matrix 
material. Here we restrict ourselves to the modification of the hardening slope alone (keeping the yield 
the same as the of the matrix) and thereby attempt to obtain the unknown artificial hardening slope E* H 
from the consideration of a simple one dimensional spring model (with each spring corresponding to a given 
suboell) which is shown in Fig. 21. 


The following discussion is limited to range 3 in Fig. 19 (i.e., the range where the effective hardening 
slope is approximate!}' constant). Assume an incremental load, AF, is applied to the spring assemblage in 
Fig. 21. This results in an incremental displacement, A u, that depends on both the geometry and material 
properties of the spring assemblage, that is, 


A F = 


. E/Af 


+ 


±±L 

E*A f 


Eh {A—Af) 


A U, 


( 2 ) 


where A is the total cross sectional area of the considered unit cell and Af the fiber cross sectional area. 
Eh and E* H are the hardening slope and artificial hardening slope, respectively. Alternatively, the global 
longitudinal stiffness behavior can be described by 


A F = 


i 


A u 


(3) 


where E € fj f is the effective longitudinal hardening slope. The combination of Ecf>. (2) and (3) , and the 
introduction of the following abbreviations, 



3 


eff _ 

r,H ~ 


" H 


Eh 



V f = ^l 
f Al 


(4) 
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A F 



Figure 21: Spring assemblage model 


lead to the final equation for the artificial hardening slope: 


E* h S i.H^ ~ M [ S rM + if *) 


(5) 


Note, the onh unknown in this equation is the effective longitudinal hardening slope E e J* which for example 
may be determined from the semi-analytical shear-lag model (.see Appendix B) or obtained from some other 

reference data. . „ 

A comparison of the FE-, GMC-. GMC-bilinear-modified-model response lustories are shown in Figs. --- 
25. Evidently, as one might expect, the CMC-model, modified by using an artificial bilinear matrix material 
defini tion in the intermediate subcell, results in a more accurate longitudinal simulation than the previous 
standard GMC-model (see Fig. 14 and that modified using an artificial elastic matrix material definition). 
With this approach, reasonably accurate results can be produced when range 2 Is small, typically for high fiber 
aspect ratios and low stiffaess ratios. Alternatively, for situations involving high stiffness ratio materials, 
low volume fractious and high aspect ratios (e.g., a,. = 25, see Fig. 24) simulations usmg the modified 
GMC approach give the largest error. This is due to the fact that these conditions induce the smallest area 
(highly loaded) subcell between fibers and thus maximizes the flow between inclusions: thereby producing 
the most compliant overall response. Clearly, only a small difference exists between FE and GMC results in 
the extensive plastic range. This Ls due to the fitting procedure utilized for the simple spring model. Note 
that, again all response histories are conservative when compared with FEA simulations. 

More detailed studies were not performed at this time using a more complex artificial matrix material 
definition, since the actual usefulness of such an ad-hoc approach is debatable due to the lack of any real 
physical meaning. Furthermore, a new higher order GMC theon (appropriately termed - high fidelity 
GMC), which provides the necessary normal and shear coupling, has recently been developed by [Aboudi et. 
al.(2001)l. This new formulation, although more computationally demanding than the standard or modified 
GMC (but significantlv faster than the FEA approach) should provide the required balance between speed 
and accuracy when attempting to simulate DRX material behavior. Consequently, further analysis using 
the GMC will be reserved until later when the theoretical extension of the present elastic continuously 
reinforced high fidelity GMC version (.see [Aboudi et. al.(2001)]) to that of discontinuous and inelastic 
material beha\ior, Ls accomplished. 
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Figure 22: Comparison of FE- and GMC-modified model ( $ r = 3, Vf = 0.2) 



Figure 23: Comparison of FE- and GMC-modified model simulations (s r =3, Vf — 0.5) 
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Figure 24: Comparison of FE- and GMC-modified model simulations (s T = 30, V f 



Figure 25: Comparison of FE- and GMC-modified model simulations (s r = 30, V) 
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5 Conclusions 


This paper has focused on evaluating the predictive ability of the generalized method of cells, in the 
context of discontinuously reinforced composite (DRX) material behavior; specifically aligned short fiber 
composites. CMC simulations were compared with: i) simulations using other micromechanical based mean 
field models and finite element (FE) unit cell modeis found in the literature for the elastic range, as well 
as ii) finite element unit cell model and a new semi-analytical elasto plastic shear lag model in the inelastic 
range. Results indicated that the GMC is full} capable of accurately analyzing the elastic behavior of DRX 
materials. However, in the case of inelastic behavior, the standard GMC approach, due to the inherent lack 
of nonnal-shear coupling, was shown to be unsuitable for predicting the longitudinal elastoplastic behavior 
of aligned short fiber reinforced composites. However, a simple modification of the GMC-model, through 
the introduction of an "artificial” bilinear matrix material definition within the subcell(s) residing between 
inclusion (fiber) ends was shown to improve the overly conservative GMC predictions within the elastoplastic 
regime significantly. This ad-hoc modification to the GMC unit cell definition would allow- reasonably 
accurate multiaxial simulations to be made. The case when both the stiffness and fiber aspect ratios are 
high being the least accurate. Finally, it was shown herein that an elastoplastic shear lag model can 
provide identical uniaxial results to that of the FE unit cell calculations. For truly consistent and accurate 
results it is recommended that a high fidelity micromechanics approach be used when examining inelastic 
discontinuous!}' reinforced composite behavior. 
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6 Appendix A: Generalized Method of Cells Formulation 

CMC can predict the elastic and inelastic thermomechanical response of both continuous and discontin- 
uous multiphased composite materials with an arbitrary internal microstructure and reinforcement shape. 
It U a continuum-based miciomechanics model that provides closed-form expressions for the macroscopic 
composite response in terms of the properties, size, shape, distribution, and response of the individual con- 
stituents or phases that make up the material 4 . These constituent materials can be represented using any 
elastic and/or inelastic deformation and life (e.g., continuum damage mechanics fatigue model) model. The 
periodic nature of composites typically allows identification of a repeating unit cell that can be used as a 
building block to construct the entire composite. The properties of this unit cell are thus representative of 
the properties of the entire assemblage once the subcell dimensions and appropriate .subcell materials are 
selected. Thus, unidirectional long-fiber composites, short-fiber composites, porous materials and laminated 
materials can all be modeled as special cases. 


6.1 Model Description 

Consider a composite material with a periodic structure whose repeating volume element consists of 
N a x N 3 x A'-, rectangular parallelepiped subcells. The volume of eadi one of the subcells is d a h 3 C ^ . where 

a, 3 and y are r unnin g indices: a = 1 N a : 3 = 1 A3 : •) = 1, .... Av,., in the x\ , 12 and 

directions, respectively. The total volume of the repeating volume element is dhC where d = X! 0 =i , 

In Fig. 26. an example is shown for a repeating volume element with N a = 2. 


h = , £-££1 

N 3 = 3 and A^ = 3. 


L 



Figure 26: Triply periodic repeating unit cell, illustrating suboell dimension nomenclature. 

By an approximate micromechanical analysis of the detailed interactions of the subcells of the repeating 
volume element, overall constitutive relations which govern the effective behavior of the multiphase elasto- 
plastic composite, can be established. This analysis relies on the requirements that static equilibrium of the 
materials in the various subcells are ensured, and that continuitv of the displacements and tractions between 
neighboring subcells within the repeating volume element, as well as between neighboring repeating volume 
elements is satisfied on an average basis. 

4 Note, each of the subcells can be comprised, in general, by an elastic-viscoplastic temperature -dependent material. 


N AS A/TM — 200 1-211165 


27 






As the average behavior of the composite Is sought, it is sufficient to consider a first order theon in which 
the displacements in the subcells are expanded linearly in terms of the distances from the center of 

each suboell. i.e., in terms of the local coordinates x[ a> , x 3 \ and xSK Thus, the first order expansion in 
the subcell (a 3 *) ) is represented as, 


(a 3*)) 




ol a3l) 


+ 2?X? a "+sW' > 


i = 1,2.3 


(A-l) 


where iv (x) are the displacement components at the center of the suboell, and o ( . a 3ri v <c ' 3 '" and 
V; are micro-variables that characterize the linear dependence of the displacement u i a3 ' l) on the local 

coordinates Xj \ x 2 \ xS ' . In Eq. (A-l) and the sequel, repeated Greek letters do not imply summation. 
Note that due to the linearity of Eq. (A-l) , static equilibrium of the material within the suboell (a 3-)) Ls 
ensured. 


The components of the small strain tensor are given by 


(a3'i) (a 5,) Jj (a3-)h . . , _ - 

( ij =2 W u j +Vj u i ) *,.7 = 1,2,3 

where = d/dx\ a) . ob = d/dx 3 ' and dz = d/dxj } . 

The volume averaged total strains and stresses in the composite are expressed, respectively, as 


(A-2) 


AL A'j 




and 


o=l 3—1 ->=1 


A„ A : 3 A', 


!j dhC T T T d “ h3C ~> n ij 




(A-3) 


(A-4) 


a = l 5=1 -) = 1 


with the relationship between the averages of the stress rr\j 3 ^\ total strain l\j 3l \ inelastic strain 

and thermal strain z- ~ a ij 1 AT ( where cv-J 3l ' Ls the coefficient of thermal expansion tensor, and AT 
is the temperature deviation T — T# from a reference temperature T^) being given by the classic Hookean 
constitutive equation 


-(a5-)) _ s-iiaSi) -I (a 3-)) -T(c\3y) , 

a ij “ U ijkC \ € kC ~ € kC ) 


(A-5) 


where C\j 3 p is the elastic stiffness tensor of the material. 

It has been shown ([Paler and Aboudi (1992)] and [Aboudi(1995)]) that by employing the displacement 
and traction continuity conditions at the interfaces between the suboeUs of the repeating volume element, 
and at the interfaces between neighboring repeating volume elements, it is passible to eliminate the micro- 
variables and obtain, via a smoothing operation (homogenization), a set of continuum equations that model 
the overall behavior of the multiphase short-fiber composite. This was achieved by establishing relationships 
which connect the microstrains at the suboelLs to the total overall macros trains in the composite via the 
appropriate concentration tensors, A and D. 


6.2 Overall Thermomechanical Constitutive Law 

Given these concentration matrices, A {Q ^ and D { ° ■■*>>, expressions for the average total strain and 
stress, respectively, in each subcell can be constructed; that is, 

e (o ^) - A ia3l) e + D {a3l) (ei + ej) 

and 

& {a 3 ~i ] = C (a3l) jA (a>>) e + D {a3l) (ei +ef) - (e I{a3 ^> + e T (^-r) 


(A-6) 

(A-7) 
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in terms of the uniform overall strain e (the applied macrostrain) and the subcell inelastic and thermal 


strains, where 


A = 


A a/ 

A c 


-.-1 


0 

J 


D = 


A m 

Ac 


-1 r 


Am 

0 


and A ; \/ involves the elastic properties C' :a3r) of the subcell material, A c the geometric dimensions ofthe 
repeating cell only, and J Ls a matrix containing the imposed average (composite) strains. Note A ° 
and D 1 " 5 "* 1 are square matrices of 6.Y a order. Similarly, notice that the same matnoes D ° are 

operating on both the inelastic and thermal strains of the subcells. 

It was shown by Dvorak and Benveniste (1992) and Dvorak (1992) that for any representative volume 
element under a uniform overall strain e and temperature change AT, which contains a piecewise uniform 
distribution of thermal and inelastic fields (eigenstrains) associated with the applied loading e . the averages 
of the local strains can be expressed in terms of the mechanical and eigenstrain concentration tensors. These 
tensors depend on the local elastic moduli, and on the shape and volume fraction of the phases and are 
therefore constant. Equation (A-6) Ls consistent vvdth this representation of the average subcell strain, where 
the constant mechanical and eigenstrain concentration tensors are given by and D ° respectively. 

Consequently, the following effective elastoplastic thermomechanical law of the composite can be estab- 

lished 

a = B*(e — c 1 — e T ) ( A - 8 ) 

where the effective elastic stiffness tensor, B*, of the composite is given by 


B- - Y. EEW 

ani a = l 3=1 -)=1 

and the composite inelastic strain tensor is defined as 

-e 1 = ~ B '~- jr - e T{a3 ^) 

dhC a=l 3^\ v=l 

and the average thermal strain tensor as 

= zEzl £ Y^jrd a h 3 c, C {a3 -"(D {a - « t(o3t) )) 

0=1 3=1 ->=1 


dhC 


(A-9) 


(A-10) 


(A-ll) 


, strains being defined as: 



f ^/(A r cvAj^)) 

(A-12) 

= (c TiUX K. 


(A-13) 


The effective coefficient of thermal expansion vector, a* of the composite, is given by 


a’ = e T /AT (A-14) 

The above micromechanically established overall thermo-inelastic constitutive law (A-S) Ls valid for any 
type of thermomechanical loading (i.e.. any combination of normal, shear and thermal loadings). A significant 
advantage of this constitutive law stems from the fact that it does not rely on any symmetry conditions that 
may exist under certain types of applied loadings. Thus in the implementation of this law, the question as 
to whether such symmetry conditions exist or not Ls irrelevant. 
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7 Appendix B: Semi-analytical Axisymmetric Elastoplastic Shear 
Lag Model 

Given the elastoplastic investigations under taken in this study, a simple shear-lag model for aligned 
short fiber composites Is introduced. As mentioned in the background section, shear lag models were the 
first micromechanics models considered for aligned short-fiber composites. Despite some serious theoretical 
flaws (e.g., the ability to predict only the longitudinal modulus £ji within the elastic range), shear lag modeLs 
have enjoyed enduring popularity; probably due to their algebraic simplicity and their physical appeal. In 
this paper, we have purposefully limited our investigation to the longitudinal elastoplastic stiffness behavior 
so that the proposed modified shear lag model can be appropriately utilize. 

Following Cox (1952), the shear lag analyses are focussed on a single fiber of length If and radius 
which is encased in a concentric cylindrical shell of radius R. Such a configuration then leads to the kev 
assumption of a shear lag model, which is that the shear stress between fiber and matrix (r ri , r.^.z. . . 
cylindrical coordinate system) Is proportional to the difference in displacement between the fiber surface and 
the outer matrix surface, i.e.. 

T r .(z) = ~ w ( r f' z )) ( B_1 ) 

where H is a oonstant and can be written as 


H = 


ln( R/r f )-V 


(B-2) 


It remains to choose the radius R of the outer matrix cylinder. Several choices have been used (see e.g. [Tucker 
and Liang (1999)]). For our purposes it appears most realistic to let R be equal to a radius corresponding 
to the width of the FE or GMC unit cell. The above equations are based on a one-dimensional analysis of 
a fiber surrounded by a matrix layer (i.e., the Poisson effects are neglected (taken to be zero)) and assumes 
linear elastic material behavior. To extend this analysis into the elastoplastic regime of aligned short fiber 
composites, the existing shear-lag model must be modified so as to: 


1 . account for elastoplastic material beha\ior and 

2. the appropriate Poisson effect. 

The elastoplastic material behavior is simply introduced by replacing G m by the plastic tangent shear 
modulus Gh (“shear hardening slope") and an incremental formulation of Eq. (B-l), i.e., 


dr ri (z) = ~-(dw(R,z)-dw(r f ,z)j 


(Bril) 


vvith 


H = 


2ttG h 


(B-4) 


In (R/r f ) - 1 

This new “shear hardening slope" can then be obtained from the P randtl-Reass-equat ion. For a pure 


shear load in the r- 2 -plane this equation can be written as 

1 3 dr r& 

dt rz = —— r/r rt + =— tv- 


2 G v 


2 r rz E H 


and leads fin all} to the shear hardening slope 

dr vz 


Gh = 


rz GmEfi 

s — = const 


de rz Eh + 3Gn 


(B-5) 


(B-6) 
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Axial-symmetry 


Matrix 



Figure 27: Sketch of the semi-analytical shear lag model 


where Gh is constant because Eh is a constant. Note, within the plastic range a fitting parameter %) is 
also introduced. This fitting parameter will account for effects w hich can not be captured by such a simple 
model. 

The Poisson effect is taken into account by considering an axis}Tnmetric model (shown in Fig. 27 using 
an exploded view') with appropriate boundary conditions on the concentric cylindrical matrix shell. 

The model posses the following characteristics: 

1. The cross section of the axisvmmetric model is divided into four regions (/, //, ///, IV). Further, 
regions III and IV are subdivided into three subregions each, e.g., Ilia , Illb, IIIc and IV a. IVb , 
IV c (see Fig. 27). 

2. Within all regions only normal stress is allowed (i.e., shear stresses are zero) such that rectangular 
cross sections of individual regions remain rectangular after loading. 

3. The longitudinal stress are constant within regions I and II and are piecewise constant within regions 
III and IV (i.e., constant within each subregion). 

4. The points 1-1\ 2-2 \ etc. are connected with non-linear truss elements which characterize the shear 
forces between the fiber and the matrix. 

5. The stiffness behavior of each non-linear truss element, i, is based on Eq. (B-l), and follows from 

F l hruss=T r .2nr f li=H(uj(R,z)-w(r f ,z))li r = C} r ^ s (w{R,z) - w(r } ,z)) (B-7) 

where, e.g., for i = 1, l\ is the distance between points 1 and 2 in Fig. 27 and represents a kind of 
shear-influence length for truss 1. C} russ represents the stiflfoess of truss 1 between points 1-1’. The 
material behavior of the truss element is bilinear and is characterized by a ‘Yield force" (which is based 
on Trescas vield condition for pure shear within the interface) and can be w ritten as 

^trusa.Y ^-2 wr f l l T (B-S) 
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as well as by a stillness hardening slope w hich can be written as 


C'trussM — *VH1 1 T . 


(B-9) 


The number of introduced truss elements depends on the desired aspect ratio. For a r — 5 three truas 
elements were used, whereas for a r = 25 six trasses were applied. 

6. Finally, in addition to the symmetry conditions shown in Fig. 27, symmetry Is also required on the top 
and right outer surface of the axis}mmetric model (see Fig. 27). 

The introduced semi -analytical shear lag model is solved numerically asing the finite element method. 

The so called semi-analytical shear-free model is obtained by assuming that the truss elements are 
characterized by zero stiffness. Such a model represents a standard GM C-model. 
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